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We use the Pade-Z2 stochastic estimator for calculating Tr log M to compute quark loop corrections to quenched 
QCD. We examine the main source of error in this technique and look at a way of controlling it. 



1. Introduction 

In lattice QCD calculation, evaluating the 
fermion determinant is computationally expen- 
sive and as a consequence det M is often not 
evaluated and is simply set equal to a constant. 
This is called the quenched approximation and 
amounts to neglecting internal fermion loops. 
The goal of this work is to study the systematic 
error introduced by quenching and to compensate 
for quenching errors by constructing an expan- 
sion in quark-loop count and using it to improve 
Monte Carlo estimators of Wilson loops measured 
on quenched gauge field configurations. We ex- 
amine a source of statistical error and a method 
for its partial alleviation. 

2. Re-vifeighting existing quenched QCD 
configurations 

With the fuU QCD action 

s[u, V, ^] = sYJ:Zn [u] + s„ [f/, V', , (1) 

the expectation value of some operator, O is given 

by 

(O) = ^ j [dU][d'iP][di,]Oe-^^^^^'^^ 

= ^ [ [df/]OdetM[i7]e-^«""«=[^] , (2) 



where 



By rewriting ^]^;'-^„[C/] as 



SYJ:^Ton[U] = nylogdetM[C/] 



= n/TrlogM[C/] 
- SS[U], 

and dividing through by e'^^^\ (||) becomes 



{O) = 



(4) 



(5) 



Expanding the exponential in the fluctuations in 
5S leads to 



+ (06S) - iO) (SS) 

\ '^S = n \ 'i5S = ^ '153=0 



(6) 



(3) 



Here, the subscript on signifies the 

quenched expectation value of O. 

3. Pade Z2 stochastic estimator for evalu- 
ating TrlogM 

In order to make use of (H) one needs to be able 
to find 5S. We use the Pade-Z2 method to find 
TrlogM as per §. 

The Pade approximant for Tr log M can be 
written as 

K 

TrlogM « boTrl + ^ 6fc • Tr(M + Cfcl)"'. (7) 
fe=i 
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Tr(M + CfcI) ^ can be estimated using noise Z2 
vectors, rj^ . 



1 ^ 



(8) 



where — (M + Cfcl) ^t;-' are the sohitions of 



(M + Cfcl)e 



k,3 



(9) 



where j = 1, • • • , L and k — I, . . . , K . The vari- 
ance of these estimators can be greatly reduced 
by subtracting suitably chosen traceless matrices 

Tr(M + ca)"' « <vH{M + Ckir' (10) 
p 

p=l 

where the Xp^k are chosen to minimise the vari- 
ance of the estimator. The Q^^-* used are 



Q 



(p) 



(1 + Ck)P+ 



j-(DP-TrDf) 



where M = 1 — kD. The odd powers of D are 
traceless, but the even powers greater than two 
are not. Tr and Tr are calculated from 
4 and 6 link loops respectively. For this work 
we expand to order K = 11, using L ~ 10 noise 
vectors and P = 8 traceless subtraction matrices. 
The variational procedure to set Xp^k was found 
not to be necessary and so we set them equal to 
unity. To solve (^) we used the MR'^ algorithm 
from |2] and the odd-even preconditioner in from 



4. Application to Wilson loops 

We apply the preceding numerical techniques 
to finding improved estimates of Wilson loops. 
Measurements were made on two sets of 10^ gauge 
field configurations; a set of 103 quenched QCD 
configurations and a set of 95 full QCD configu- 
rations, generated using the MILC collaboration 
software 

Both have /3 = 5.436, and the fuU QCD con- 
figurations have K ~ 0.16 and nj — 2. Fig. |l| 
shows the results for various Wilson loops with 
corrections done to first order as per (||) In all of 
Figs. |l] ~ ^ the loop number on the x~axis labels 
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Figure 1 . Effect of correcting Wilson loops by ex- 
panding to first order in 5S[U] — n/TrlogM[C/] 
as per (^. 



Wij from left to right in the following order: ^ , 

wwwwwwwww 

^1 10 ' ^2 5 ■ Figs. |l|, H and ^, points have been 
offset slightly for clarity. 



5. A method for controlling the statistical 
error 

A method for controlling the statistical error 
in (^) is suggested by looking at the form of the 
correlator between the diagonal elements of log M 
and the operator, O, being measured. Given that 
O and 5S can be written 



= ^0, , 5S^y^{5S) 



y 



y 



(11) 



then the first order correction to O can be writ- 
ten as 



HO) 
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\ X I \ y I 

J2{0.{SSy-{SSy))) 



(12) 



x,y 



If we write /(r) 
then 



OxSSy), where r = \x ■ 



Ox5Sr 



drdflr^fir) 



(13) 



E 

Fig. H shows /(r) //(O) and r^f{r) //(O) as a func- 
tion of for the case O = W^^ . The former shows 
that the two operators are correlated, as would be 
expected. The latter shows that most of the sig- 
nal for 80 is at small r and but most of the noise 
is at large r, and that for large lattices this has 
the potential to swamp the signal. One way to 
address the problem is to cut the integral off at 
some rm.aT.- 



{60)= J2 (^-^'^; 



(14) 



This has been done using rmax = 4, and the re- 
sults are shown in Fig. ^. The cut off scheme 
reduces the size of the statistical error but also 
looses some signal. 

The Monte Carlo estimator for (O) is 

1 ^ 

«'^^))^1^E E ^"-^^y- (15) 



n— 1 \^ — v\ 



An efficient way to do this is with Fourier con- 
volution. Given a windowing function W{x — y), 
where 



W{x -y) = 



1 ii \x — y\ < r 



otherwise 
then (|l^) becomes 



(16) 



n—l x,y 

- ^EE^" (17) 

n— 1 X 

X [T-'[T{SS-}HW{x,y)}}]^ 



where J- and J-~^ are defined in the usual way 
and implemented as fast Fourier transforms. Us- 
ing this method to compute (SO) also makes it 
easy to choose a more elaborate windowing func- 
tion should this be desired. 
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Figure 2. Contribution to {60) of correlator be- 
tween iSS[U])y = (n/ \ogM[U])yy Bud {W,,)a: to 

m 



Approximating full QCD by shifting the 
gauge coupling and expanding Tr log M 



Following ||5| and identifying 

SS = 71/TrlogM - SpSgauge , 



(18) 



we again look at corrections to Wilson loops. 
These are measured on a set of 200 10'' full QCD 
configurations with /? — 5.679, k = 0.16 and 
Uf — 2 and also the set of quenched QCD con- 
figurations used earlier. The results are shown 
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Figure 3. Corrections to Wilson loops relative to 
the full QCD value comparing the effectiveness 
with and without the cut off scheme. 



measurements. 
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in Fig. ^ and are in accordance with ||]. A no- 
table difference between the two is the computa- 
tional savings afforded by the choice of algorithm 
for Tr log M. Use of unbiased subtracters means 
that only 10 noise vectors are needed. In addition 
to this we do not have to fix the configurations to 
Landau gauge. 

The previously mentioned method of cutting 
the integral off to reduce the statistical error turns 
out not to work for 6S defined in (|l8|). The con- 
tributions from Tr log M and 5j3Sgauge are nearly 
equal but are opposite in sign and so the total 
correction is small as is shown in Fig. |^ 

Introducing a cutoff reduces the contributions 
from Tr log M and 5 (5 S gauge by different amounts 
which produces a large variation in the total cor- 
rection. This is shown in Fig. |^. 

We have also looked at corrections to two point 
functions but the noise is such that this method 
cannot be used as is to provide improved mass 



Figure 4. Relative corrections to Wilson loops by 
expanding to first order in SS = TiyTrlogM — 
SpSgauge as per (|l|) 



7. Conclusion 

A method of improving estimates of observ- 
ables measured on quenched QCD configurations 
was tested on some sets of small (10*) config- 
urations. Expanding to first order in 5S = 
n/TrlogM gave moderate improvement to esti- 
mators for Wilson loops. Analysis of the source of 
the statistical error hinted at problems that the 
method would have on larger lattices and with 
two point functions, and also suggested a tech- 
nique for its partial alleviation. The technique 
does have the potential to address these prob- 
lems but at the price of reducing the size of the 
correction. 
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Figure 5. Relative corrections to Wilson loops 
from each component of 6S both separately and 
in combination. 



Figure 6. Relative corrections to Wilson loops 
from each component of SS using the cut-off 
scheme to control the statistical error, and with 
these components combined. 



It was shown in |^ that combining the above 
method with an appropriate shift in the gauge 
coupling can provide more accurate corrections to 
Wilson loop values. We obtained the same results 
using a more economical method of calculating 
TrlogM. 
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